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ABSTRACT 



CNJ ' We discuss the nonlinear development of the isobaric mode of thermal instability (TI) in the context of 

^ , the atomic interstellar medium (ISM), both in isolation and in the presence of either density or velocity 

^ ' fluctuations, in order to assess the ability of TI to establish a well-segregated multi-phase structure in 

I the turbulent ISM. The key parameter is the ratio of the cooling time to the dynamical crossing time 77. 

ly-^ , First, we discuss the degree to which the condensation process of large-scale perturbations generates large 

■ velocities, and the times required for them to subside. Using high-resolution simulations in ID and fits to 
\ recently published cooling rates, we find that density perturbations of sizes > 15 pc in media with mean 

■ density ~ 1 cm~'^ develop inflow motions with Mach numbers larger than 0.5 and a shock that propagates 
^ ' outwards from the condensation, bringing the surrounding medium out of thermal equilibrium. The time 

, for the dynamical transient state to subside ranges from 4 to 30 Myr for initial density perturbations of 

^ ■ 20% and" sizes 3 to 75 pc. By the time the condensations have formed, a substantial fraction of the mass 

' is still traversing the unstable range. Smaller (0.3 - 3 pc) perturbations may condense less dynamically, 

^ID . and reach nearly static configurations in shorter times (e.g., ~ 3.5 Myr for perturbations of ^ 0.3 pc), but 

' they may be stable if they have a turbulent origin (see below). We thus suggest that, even if TI were the 

\ sole cloud- forming agent in the ISM, clouds formed by it should be bounded by accreting gas traversing 

■"""^ . the unstable range, rather than by sharp transitions to the stable warm phase. Second, we discuss the 

' competition between a spectrum of density perturbations of various sizes. We empirically find that, in 

order for small-scale perturbations not to alter significantly the global evolution, progressively larger 
Q ■ values of -q are necessary as the initial spectrum becomes shallower. Finally, we discuss the development 

I of the instability in the presence of random velocity forcing, which we argue is the most realistic way to 

C/3 . emulate density fluctuation production in the actual ISM. Such fluctuations are quasi-adiabatic rather 

' than quasi-isobaric in the large-ry limit, and trigger the wave mode of TI, rather than the condensation 

^ \ mode, being stable to first order. Indeed, we find that the condensation process can he suppressed for 

arbitrarily long times if the forcing causes a moderate rms Mach number (> 0.3) and extends to small 
enough scales or occurs in low enough density environments that the turbulent crossing time becomes 
smaller than the cooling time at those scales. We suggest that this mechanism, and the long times 
required to evacuate the unstable phase, may be at the origin of the relatively large amounts of gas mass 
in the unstable regime found in both observations and simulations of the ISM. The gas with unstable 
temperatures is expected to be out of thermal equilibrium, suggesting that it can be observationally 
distinguished by simultaneously measuring two of its thermodynamic variables. We remark that in the 
(stable) warm diffuse medium 77 is large enough that the response to velocity perturbations of scales 
up to several parsecs is close to adiabatic, implying that it is relatively weakly compressible, and thus 
consistent with recent observations that suggest a nearly Kolmogorov power spectrum in this medium. 

Subject headings: instabilities — ISM: clouds — ISM: kinematics and dynamics 



1. INTRODUCTION 

The fact that the neutral atomic interstellar medium 
(ISM) is most likely thermally bistable (Zeldovich & 
Pikel'ner 1969; Field, Goldsmith & Habing 1969; Wolfire 
et al. 1995) has had a great impact on our picture of the 
formation of interstellar structure. In particular, the un- 
stable cooling process has been thought to be responsible 
for the collection of mass into HI clouds of sizes ~ 1 pc, 
densities ~ 10-100 cm^^, and temperatures ^ 100 K (e.g.. 
Goldsmith 1970; Schwarz, McCray & Stein 1972; Parra- 
vano 1987; Lioure & Chieze 1990; Koyama & Inutsuka 



2002). In the classical model of McKee & Ostriker (1977), 
the cold and warm neutral media are two distinct phases, 
with clouds forming in dense shells created by shock com- 
pressions, aided by the thermal instability (TI hereafter), 
and later confined by the warm medium's thermal pres- 
sure. The triggering of TI in the warm medium by strong 
compressions and the subsequent segregation in two sep- 
arate phases has been recently studied by Hennebelle & 
Perault (1999, 2000) and by Koyama & Inutsuka (2000). 

However, over the years, a more complex picture of the 
ISM has emerged, driven in part by the advent of realistic 
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numerical simulations of the ISM that have highlighted the 
role of inertia and advection (i.e., transport) by the veloc- 
ity field. In this context, the role of the thermal instability 
in the overall dynamics of the ISM is still uncertain. Simu- 
lations of the turbulent ISM driven by supernovae or stel- 
lar winds suggest that the complex structure of the ISM 
as well as HI clouds can be the result of random gas mo- 
tions, at least in regions with densities below 10-100 cm~^, 
with no need for thermal instabilities (e.g. Bania & Lyon 
1980: Rosen & Brcgman 1995; Passot, Vazqucz-Semadcni 
& Pouquet 1995; Korpi et al.^ 1999; see also Vazquez- 
Semadeni & Passot 1999 and Vazquez-Semadeni 2002 for 
reviews; but see Wada, Spaans & Kim 2000; Koyama & 
Inutsuka 2002; Kritsuk & Norman 2002 for the alternative 
view that TI powers ISM turbulence). 

Within this scenario, Vazquez-Semadeni, Gazol & Scalo 
(2000, hereafter Paper I) studied the probability density 
distribution (PDF) of mass density in low resolution (128^) 
2D simulations of the thermally bistable ISM at large 
scales (1-kpc box), finding that the PDF may not reflect 
the thermally unstable nature of the gas in the presence of 
magnetic fields. Galactic shear, the Coriolis force, and tur- 
bulence driven by stellar-like sources located at the density 
maxima (clouds). Moreover, simulations at higher resolu- 
tion (512^) with both the same box size (1 kpc; Gazol et 
al. 2001, hereafter Paper II), and smaller ones (150 and 
10 pc; Vazquez-Semadeni et al. 2002, hereafter Paper III) 
have temperature histograms that indicate the existence 
of significant amounts (~ 50%) of gas with temperatures 
in the unstable range. The same result could be seen in 
the simulations by Korpi et al. (1999), although they did 
not discuss it, and has later also been found in simula- 
tions of the late nonlinear stages of the pure instability at 
smaller scales (6-pc box, at 256^ resolution) by Kritsuk & 
Norman (2002). This is in agreement with observations 
by Dickey, Salpeter & Terzian (1978), Kalberla, Schwartz, 
& Goss (1985), Spitzer & Fitzpatrick (1995), Fitzpatrick 
& Spitzer (1997), and Heiles (2001), as well as previous 
theoretical suggestions (e.g., Lioure & Chieze 1990; Nor- 
man fc Ferrara 1996). Thus, the study of TI in the context 
of fluctuating media becomes important, in particular to 
understand the origin of the large fractions of gas mass in 
the unstable regime. 

At the large scales studied by Paper I, the condensation^ 
process leading to the formation of clouds and intercloud 
"voids" becomes highly dynamic (i.e., develops large ve- 
locities in the nonlinear stage) . The parameter controlling 
the degree to which the condensation process is dynamic 
is the ratio (ry) of the cooling to the representative (sound 
or bulk) dynamical time. Although it is customary in the 
literature to discuss TI in terms of the perturbation wave 
lengths (or wavenumbers; e.g. Field 1965), in this paper 
we prefer a description in terms of r] because it allows the 
consideration of perturbations of a given physical size (say, 
in parsecs), but that may correspond to different values of 
T] because of, say, they occur at different mean densities 
(i.e., equilibrium cooling rates). 

It is well known (e.g.. Field 1965; Murray & Lin 1989; 

^ The simulations of Korpi et al. actually use a cooling function that imply a thermally unstable range between the cold and warm gas if no 
heating is present (see, e.g., Murray Sz Lin 1989; Burkert & Lin 2000). Although no global heating was used in those simulations, the presence 
of local stellar heating causes the gas to be thermally stable while subject to heating, and unstable while cooling. 

^ For consistency with the nomenclature in the literature, we maintain the name "condensation" even for highly dynamic cases in which 
"accretion" or "collapse" might be more appropriate. 



Burkert & Lin 2000) that perturbations with rj > 1 (cor- 
responding to small-scale perturbations at a fixed cooling 
rate or low cooling rate at a fixed length scale) evolve al- 
most isobarically, and thus quasi-statically. Instead, per- 
turbations with T] < 1 develop large pressure gradients, 
and thus evolve highly dynamically even under the isobaric 
mode of the instability, although they take longer times to 
grow (e.g.. Field 1965; Balbus 1995; Meerson 1996). Thus, 
Paper I suggested that large clouds formed by TI should 
be bounded by accretion shocks (rather than by contact 
discontinuities), in a similar manner to that discussed by 
David, Bregman & Seab (1988) for TI in lO^-lO^-K gas, or 
by Rees (1977) in the case of gravitational contraction. It 
is then important to determine whether large, static, ther- 
mal pressure-confined clouds ever develop and, if yes, how 
long does it take for the dynamic situation to fade away. If 
this time is long compared to estimated cloud lifetimes [see 
Blitz & Shu (1980) and Ballesteros-Paredes, Hartmann & 
Vazquez-Semadeni (1999) for two different possibilities in 
the case of molecular clouds], or to the time scales for 
major disturbances to pass through the cloud (e.g., Korn- 
reich & Scalo 2000), then the static configuration should 
not be expected to occur in general in the ISM. Unfortu- 
nately, Paper I did not proceed to give a definitive solution 
to this question because of limitations of the numerical 
scheme used there. Most previous discussions of the time 
scales to reach a developed cloud-intercloud medium have 
been based on the assimiption of isobaric evolution (e.g., 
Parravano 1987), which, however, is strictly valid only in 
the limit of r] » 1. Goldsmith (1970) did consider the 
limit 77 -C 1, and mentioned that very long times are re- 
quired for the quasi-static regime to be reached, such that 
his simulations (at very low resolutions) never reached it 
(see also Schwarz et al. 1972). 

Small-scale perturbations, on the other hand, are good 
candidates for rapid, yet not highly dynamic condensa- 
tion, and constitute the paradigm of cloud formation via 
the isobaric mode of TI. This brings up the questions of a) 
whether, in the presence of a spectrum of density fluctua- 
tions, the small-scale perturbations can use up most of the 
available mass before the large-scale ones grow, and thus 
impede their development, and b) whether the paradigm is 
maintained when the perturbations are produced by tur- 
bulent velocity fluctuations. In the regime of Paper I, in 
which the cooling times are short, turbulence forced at 
large scales cannot prevent the development of TI, and in 
fact actually enhances it, although stellar-like heating at 
small spatial scales was shown to erase the signature of 
TI in the density histogram. The latter result does not 
strictly constitute an inhibition of TI, because the stars 
form after the condensations, which in turn formed by a 
combination of the turbulent compressions and TI. Such 
condensations, in the absence of the stellar heating, would 
pile up in the thermally-stable density regime. However, 
true avoidance of the condensation may occur for random 
velocity fluctuations at the small scales because in this 
case the compressions tend to heat the gas; i.e., the per- 
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turbations are quasi-adiabatic, and can re-expand before 
they have time to cool. 

The discussion above is related to the question of 
whether the classical picture of cold, dense clouds in pres- 
sure equilibrium with their warm, dilute surroundings is a 
realistic model of the ISM, even if TI were the sole cloud- 
forming agent in the ISM (which it is not). Indeed, the 
equilibrium configuration is a possible solution of the hy- 
drodynamic equations in the presence of Tl-inducing cool- 
ing functions; it is very often used as an initial condition 
of both analytical and numerical studies, and even as a 
textbook case. However, the question is whether it is also 
a self-consistent final state of previous evolution starting 
from small-amplitude perturbations of thermally unsta- 
ble gas. Of course, the pressure equilibrium configuration 
is valid as an initial condition for numerical simulations 
or as a case study for theoretical analyses, just as Gaus- 
sian random fluctuations arc in studies of turbulence or 
Cosmology, but they are not to be interpreted as realis- 
tic representations of the ISM. For example, shock-cloud 
collision simulations (e.g., Klein, McKee & Woods 1995) 
often start with a cloud in pressure equilibrium with its 
environment, but end up with fully developed turbiilence, 
suggesting that the latter is a more likely regime for the 
ISM than their idealized initial conditions. 

In this paper we explore the nonlinear development of 
TI in the isobaric mode, in the context of the atomic ISM, 
focusing on the issues mentioned above, using high reso- 
lution numerical simulations in one dimension (ID) and 
a numerical scheme better suited for this type of prob- 
lems than that used in Paper I. We should emphasize that 
in this paper we choose to isolate the effects of the in- 
stability, and concentrate on a somewhat idealized case 
(one-dimensional geometry, simple parameterized cooling 
functions neglecting non-equilibrium chemical evolution, 
neglect of other physical agents such as the magnetic field, 
differential disk rotation, self-gravity and realistic, stel- 
lar energy injection), not necessarily representative of the 
full situation in the ISM. We make no attempt to discuss 
the fate of the clouds formed by TI in the real ISM, as 
this must be accomplished by means of simulations of the 
fully turbulent ISM in the presence of the agents men- 
tioned above (see Papers I, II and III), and is likely to 
involve dynamical instabilities whose study requires two- 
or three-dimensional simulations (see, e.g., Murray et al. 
1993; Vietri, Ferrara & Miniati 1997). Instead, our main 
goals here are simply to illustrate, even in this simplified 
case, the following features of the condensation process: 
a) the dynamical nature of the condensation process, pro- 
viding quantitative estimates of the perturbation sizes in 
the atomic ISM required for such dynamical evolution to 
occur, and of the times required for it to subside; b) the 
existence of accreting gas, which remains at densities in 
the imstable range, even after the condensation process 
has reached a slowly evolving stage; c) the competition 
between growing modes of different sizes; d) the evolution 
under the isobaric mode of TI of perturbations caused by 
random velocity fluctuations, and e) their implications on 
standard scenarios of a multiphase ISM in which most of 
the gas must reside in one of the stable phases. 

The plan of the paper is as follows. We first present 
the numerical method and discuss the resolution require- 



ments (§2). Next, we present a qualitative discussion of 
the timescales involved and the regimes that can appear 
depending on the ratio r] (§3), to then consider the evolu- 
tion of Gaussian density perturbations, the size scales at 
which shocks form at advanced stages of the "condensa- 
tion" process, and the times required to reach quasi-static 
states (§4.1). We then discuss the competition between the 
growth of small- and large-scale perturbations by means of 
simulations with an initial spectrum of density fluctuations 
(§4.2), and the evolution of turbulent fluctuations in the 
presence of TI by means of simulations including random 
forcing (§4.3). In §5 we compare our results with previous 
work, and discuss the limitations of the one-dimensional 
simulations we have used as well as some implications of 
our results on the classical picture of thermal-pressure con- 
finement of clouds, the existence of gas in the unstable 
regime, and the compressibility of the warm gas in the 
ISM. Finally, we summarize our results in §6. 

2. THE MODEL 

2.1. Numerical method and relevant parameters 

We solve the hydrodynamic equations for the logarith- 
mic density Inp, the velocity v, and the entropy per unit 
mass s. In terms of these variables, the mass, momentum 
and energy conservation equations used to calculate the 
complete evolution of a thermally unstable medium are 
D\np 



Dt 



= -V • V, 



Dv _ VP 

'Dt ~ y 

Ms 



f{r, t) + viscous force, 



T- 



Dt 



r — pA -|- viscous heating. 



(1) 
(2) 



(3) 



where D/Dt = d/dt -|- w • V is the advective derivative, 
/ is a forcing function to be specified that will be used 
only in the case of forced simulations (§4.3), and F and 
A arc the heating and cooling functions, respectively, de- 
pending on p and T. The equation of state used is that 
of ideal gases P = (7 — l)c^pT, where 7 = 5/3 is the 
heat capacity ratio Cp/c^ for monoatomic gases. The spe- 
cific entropy is related to the density and pressure by 
s = soo + Cv [ln(P/Poo) - 7ln(/9/poo)], where pqq, Pqo and 
soo are reference values for the density, pressure and en- 
tropy. 

We use simplified forms for the cooling and heating func- 
tions, consisting of a piecewise power-law fit to the "stan- 
dard" equilibrium P-p curve of Wolfire et al. (1995, here- 
after WHMTB) and a constant background heating rate 
per unit mass. This choice is at an intermediate level in a 
hierarchy of increasing complexity in the treatment of the 
thermal processes in a gas dynamical problem, that goes 
from assuming isothermal (e.g., studies of molecular clouds 
by Gammie & Ostriker 1996; Mac Low et al. 1998; Padoan 
& Nordlund 1999; Ostriker, Stone & Gammie 2001) or 
polytropic flows (e.g., Vazquez-Semadeni, Passot & Pou- 
quet 1996), to solving the heat equation using simple fits 
to the cooling and heating functions (e.g., Chiang & Breg- 
man 1988; Rosen & Bregman 1995; Vazquez-Semadeni, 
Passot & Pouquet 1995, 1996; Paper I; Burkert & Lin 
2000; Koyama & Inutsuka 2002; this paper), to follow- 
ing the detailed non-equilibrium chemistry and the cool- 
ing rates (e.g., David et al. 1988; Murray & Lin 1989, 
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1990; Kang et al. 1990; Koyama & Inutsuka 2000). For 
the purpose of studying the dynamical development of TI 
from the equilibrium state in the neutral ISM, it is not 
clear to what extent it is imperative to follow the detailed 
non-equilibrium chemistry and the resulting rates (sec §5.1 
for further discussion). Thus, the present work should be 
taken mainly as an illustration of the highly dynamical 
and complex regimes that can already appear within this 
approximation, and the caveat of the simplified cooling 
functions wc use should be kept in mind. 

Our fit to WHMTB's standard cooling curve has the 
form 

A = a,i+iT'5*.*+i for Ti<T<Ti+i. (4) 

As mentioned above, we parameterize the heating func- 
tion by r = ro(p//9o)") with a = 0, since the dominant 
heating process (photoelectric heating from small grains 
and PAHs) in the density range of interest depends only 
weakly on density, roughly as p"-^ (WHMTB). 

Wc find the following set of /3- values and transition tem- 
peratures: /?i2 = 2.12, /323 = 1.0, P34 = 0.56, 1345 = 3.67, 
and Ti = 10 K, T2 = 141 K, T3 = 313 K, T4 = 6102 
K, Ts = 10'"^ K. Note that for T3 < T < T4 the gas 
is thermally unstable under the isobaric criterion, and 
marginally stable for T2 < T < T3. For a constant heating 
of To = 0.015 erg s~^g~^, the values of the coefficients of 
the cooling functions in units of erg s~^g~^cm^K~^*'*+i are 
Ci,2 = 3.42 X 10^6^ (^2,3 = 9.10 X 10^^, C3A = l.H x 10^° 
and €4^5 = 2.00 x 10*^. The fitted curve is displayed in 
Figure 1 (c.f. fig. 3 of WHMTB). 

We will denote by Peq{p) the pressure corresponding 
to thermal equilibrium at density p, by Pmax the equi- 
librium pressure corresponding to the high-temperature 
boundary of the thermally-unstable range, and by Pmin 
the low-pressure plateau in the thermal equilibrium curve 
(see Fig. 1). As "standard" unperturbed initial conditions 
we take = 1 cm~^, representative of the mean density 
in the ISM. and T = 2400 K, the thermal equilibrium tem- 
perature corresponding to po- This state falls in the unsta- 
ble range. Other workers have chosen to apply nonlinear 
perturbations (strong compressions or shocks; e.g., Hen- 
nebelle & Perault 1999; Koyama & Inutsuka 2000, 2002) 
to linearly stable gas in order to determine whether TI 
can be triggered in such regimes by strong perturbations. 
However, since in this paper we are interested in the con- 
ditions necessary for the development of TI to be highly 
dynamic and in the competition between modes, we prefer 
to take the initial conditions described above, which, we 
remark, are close to the mean density of the ISM (see, e.g., 
Ferriere 2001). 

Wc use a sixth-order finite difference scheme for the spa- 
tial derivatives on a one-dimensional, non-uniform Carte- 
sian grid, and a third-order time stepping scheme, as de- 
scribed in Sanchez-Salcedo & Brandenburg (2001). The 
code includes a shock-capturing viscosity proportional to 
— V • V wherever V • f < 0, and zero otherwise. Shocks 
are resolved over 3-4 zones. To avoid numerical insta- 
bility due to amplitude errors at high-wavenumber modes 
of density and entropy, diffusion terms in the continuity 
and entropy equations have been introduced. Convergence 
tests (§4.1.1) have been performed to ensure that mass and 
thermal diffusion do not affect the growth rate of forming 
density peaks. In particular, the effective mass diffusion 



scale implies that below scales 10 pixels the growth rate 
starts to decrease monotonically (see also Paper III). 

Most of the simulations presented here assume a plane- 
parallel geometry. This is motivated by the need to achieve 
very high resolutions, but is justified by evidence that sug- 
gests the presence of sheet-like and filamentary morpholo- 
gies in the cold neutral medium (e.g. Kim et al. 1998). 
Although the origin of this morphology may not be caused 
by TI, we consider this geometry to be more appropriate 
than, for example, a spherical one for studying the devel- 
opment of TI in the neutral medium of the Galaxy. We 
discuss this issue further in §5.2. 

The non-uniform grid allows for enough resolution as to 
resolve the structure of the cloud and the accretion fronts. 
It also allows us to move the boundaries far away, so no 
perturbation has time to reach them. For simulations of 
forced flows, we use a uniform mesh with periodic bound- 
ary conditions. 

2.2. Resolution requirements 

One condition on the the minimum resolution necessary 

in the simulations is obtained by requiring that the final 
condensations be well resolved. Consider a domain of size 
I per dimension covered with a uniform mesh containing 
Ux grid points, in the x-direction, separated by a distance 
6x, and an initial density perturbation with a characteris- 
tic size Iq = iJjI. In ID, the final size of the condensation, 
li, will be h ~ IqPq/pi = iplpo/pi, where po and pi are 
respectively the mean and final densities. If we require 
a minimum number of zones, n, inside li, the total num- 
ber of grid points must be larger than npi/{'ippQ). For 
Pi/ Pa = 50, n = 10, and ip = 0.5, we need Ux > 1000. 
If this condition is not met, our tests indicate that a sta- 
tionary regime is reached in which mass diffusion balances 
accretion, leading to a very unreahstic situation. 

If one wants to resolve the out-of-equilibrium thermal 
structure within accretion fronts, a second condition must 
be satisfied, namely that 6x must be smaller than the cool- 
ing length. The cooling length is given by = TcooiVh, 
where TcooI is the cooling time and the front velocity. 
The condition reads > I / {Tcooi'Vfi) ■ This is a highly 
variable condition, as both TcooI and Vf^ decrease with in- 
creasing density, so the condition is most stringent at high 
densities. For TcooI ~ 0.2 Myr and Vfr = 0.1 km we 
need n,,. > 10"* for I = 100 pc and n,. > 10'"' for Z = 1 
kpc. The need for huge resolution can be alleviated by us- 
ing a non-uniform grid, with a maximum density of points 
where the condensation is expected to form. 

An alternative approach is to reduce the heating and 
cooling coefficients by some factor. This modification of 
the cooling rates is justified only if 77 -C 1 at all scales of 
interest, and at all times, since, in effect, this is equivalent 
to reducing the simulation box size. However, in regions 
where the gas is evacuated, the density decreases and the 
cooling times increases. Hence, it is not true that the con- 
dition <C 1 holds at all times in the warm diffuse gas, and 
we cannot reduce the heating and cooling coefficients by 
the same factor everywhere, but instead we would have to 
use a density-dependent function which reduces the cool- 
ing rates only at very high densities. For these reasons, 
we prefer to stick to non-uniform grids in simulations of 
single condensations, and to very-high-resolution uniform 
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grids in multi-fluctuation simulations. Note, however, that 
these high resolution requirements arise only in order to 
resolve such small-scale structures as the cloud formed or 
the thermal structure within radiative shocks. Nonethe- 
less, much lower resolutions still suffice to describe the 
global dynamics. 

3. TIME SCALES AND GENERAL FEATURES 

In this Section, we describe some basic concepts use- 
ful for understanding the remainder of the paper, includ- 
ing the characteristic time scales involved, as well as some 
features of the growth of fluctuations. We also define the 
relevant time scales explicitly in the cases of density and 
velocity fluctuations. 

3.1. Modes of the instability and the dynamical evolution 
of fluctuations 

As mentioned in the Introduction, in this paper we are 

concerned with cooling functions which imply instability 
only for the isobaric mode. In the presence of constant 
heating {a = 0), this is satisfied for /3 < 1, in which /3 
is the exponent of the temperature in the cooling func- 
tion (see §2). For our fit to the standard cooling function 
of WHMTB, this occurs in the temperature range 300 K 
< T < 6100 K, where /3 = 0.52, i.e., between the so called 
"cold" and "warm" phases of the ISM. In the presence 
of constant heating, the isochoric mode of TI is unstable 
for /3 < 0, and thus for temperatures between ~ 10^ and 
~ 10® K for classical cooling functions (e.g., Raymond, 
Cox & Smith 1976). Such temperatures are not reached by 
our simulations, since we do not include supernova heat- 
ing. Finally, it can be easily shown that in our case the 
adiabatic criterion requires /3 < 1/(1—7) —3/2 for insta- 
bility, and thus our simulations arc stable under the latter 
two criteria in the vicinity of thermal equilibrium. Note 
that, however, far from thermal equilibrium, the system 
may cool freely, and then the isochoric mode may also be- 
come unstable, as in the cooling case (i.e., in the absence 
of background heating) , in which the instability criteria for 
the isobaric and isochoric modes are, respectively, (3 < 2 
and /? < 1. 

3.2. Evolution of isobaric density perturbations 

Here we consider the growth of initial states with density 
perturbations of a given physical size, no fluid motions and 
constant pressure. Isobaric density perturbations consti- 
tute a simple case of entropy perturbations, as they imply 
a perturbation of the ratio P/p'', and can be produced 
in practice by either locally varying the cooling or heat- 
ing rates, or by velocity fluctuations well into the subsonic 
regime (sec §§3.3 and 4.3). 

Let us consider, for simplicity, the condensation process 
of a perturbation pi = p — po > 0. For a linear pertur- 
bation of wavelength A there are two relevant timescales 
at t = 0: the soimd crossing time Ts(0) = A/cs(0), 
where Cs(0) is the adiabatic sound speed, and the cool- 
ing time Tcooi (0) = Cv/ \{dpA/dT)p\, which is of the order 
of w cq/pqA, where eg is the initial internal energy per unit 
mass. The initial time scale ratio is then rjo = Tcooi/ts- 

It is well known (Field 1965; see also Shu 1992; Balbus 
1995; Meerson 1996; Burkert & Lin 2000) that entropy 
perturbations with 77 > 1 evolve almost quasi-statically. 



maintaining quasi-isobaricity throughout their evolution. 
On the other hand, entropy perturbations with 77 < 1 cool 
before they can equate pressures and therefore develop 
large pressure gradients, which then cause them to evolve 
highly dynamically even if only the isobaric instability cri- 
terion is satisfied. It is also well known that perturbations 
with T] > 1 grow on time scales ~ TcooI, because in this 
case the perturbation can only evolve as fast as it cools, 
and thus the growth rate saturates at ~ 1 /tcooI in the limit 
of small perturbation scales. Instead, when only the iso- 
baric criterion is satisfied, perturbations with 77 < 1 grow 
on a time ~ Tg, because they can only evolve as fast as 
the pressure gradient pushes the gas surrounding the con- 
densation to equate pressures (see Field 1965 and Meerson 
1996 for more details). This implies that the growth rate 
approaches zero in the limit of long wavelengths. In order 
for finite growth rates to exist in this limit, it is necessary 
to also satisfy the isochoric criterion. 

If the evolution were strictly isobaric, a condensation 
(cloud) would end up acquiring the density denoted pisob in 
Figure 1. Adopting the cooling function described in §2.1 
and a mean interstellar density of 1 cm~^ (Teq ~ 2400 K), 
Pisob = 40 cm~^. Furthermore, the condition Tg ^ TcooI 
occurs at a size scale A ~ 10 pc (hereafter denoted by 
Zeq)- Therefore, isobaric evolution cannot occur at scales 
larger than a few parsecs. Let us thus consider a pertur- 
bation (of large enough size scale) such that <C 1. Since 
Tcooi oc p~^l^ and Ts cannot decrease faster than p~^l^ 
during the condensation process, where D is the dimen- 
sion of the contraction, we conclude that <C 1 in the 
core of the cloud (i.e. at the vicinity of x ^ 0) at all subse- 
quent times, even as it shrinks. This is in sharp contrast 
with the fact that smaller scales have larger values of 77, 
under the mean conditions of the unperturbed medium. 

Relatively high local Mach numbers can be attained in 
the condensation process imder conditions with 770 <C 1. 
An estimate of the maximum velocities that can be reached 
due to pure TI can be made from Bernoulli-like arguments 
as v^^^y w 2(P„iax — -PmiiO/po, leading to fluid velocities 
Winst 4 km s"\ for po = 1 cm"-\ Pm&yi/k = 3000 K 
cm^"* and Pmin/k = 1000 K cm^'^. This velocity is clearly 
supersonic in a medium at the temperature corresponding 
to Pmin, at which Cg » 2 km s~^. With that typical veloc- 
ity, a region of size L (Lq at time t = 0) condenses on a 
characteristic timescale 

= = 1.2 f-i:^) 10® yr. (5) 

The characteristic time given by Eq. (5) can be com- 
pared to the free-fall time of a uniform cloud due to self- 
gravity, 

Therefore, self-gravity must be included either at scales 
larger than ~ 0.4 kpc or at late stages due to the density 
enhancement. In a plane-parallel condensation Tf[{t) oc 
L^^"^, so that the ratio Tg/vmst increases during the con- 
densation process. Thus, for the fiducial values we have 
used, the flow may reach the thermally stable^ regime 
and re-approach pressure equilibrium, leaving only a slow- 
accretion regime (see §4.1), without self-gravity ever be- 
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coming dynamically important. Consequently, overshoot- 
ing of density will be unable to prompt gravitational in- 
stabilities in ID experiments, but it may be important in 
more than ID (e.g., Schwarz et al. 1972; Kang et al. 2000). 

3.3. Fluctuations in the velocity field 

If, in contrast with the cases described above, an initial 
state of constant density is perturbed with velocity fluctu- 
ations (of amplitude vq), the evolution may be quite differ- 
ent. The generalization of the instability analysis to a non- 
static initial state was performed by Hunter (1970, 1971). 
Velocity perturbations as triggers of TI in initially sta- 
ble media have been investigated by a number of authors 
(e.g., Murray & Lin 1989; Kang et al. 1990; Henncbelle & 
Perault 1999, 2000; Koyama & Inutsuka 2000), but here 
we are interested in the opposite effect, i.e., whether per- 
turbations induced by velocity fluctuations in initially un- 
stable media can actually be stable, as first suggested, but 
not shown, in Paper I. A thorough discussion has been 
presented in Paper III, and here we just repeat the main 
points. 

In the presence of velocity fluctuations, a new charac- 
teristic time scale appears in the system, namely the bulk- 
velocity crossing time, Tu = X/u, where u is the character- 
istic velocity difference across scale A. Thus, the dynamical 
time Tdyn to use in r] should be chosen as the minimum of 
the sound and the bulk velocity crossing times. 

Compressive motions in the absence of radiative cooling 
cause adiabatic perturbations in which the density and 
pressure fluctuations have the same sign, because of the 
PdV work done on the affected fluid parcel. In the pres- 
ence of cooling, the actual behavior depends on (the re- 
deflned) rj. If the cooling time is much longer than the 
characteristic time of the compression (i.e., if r/ ^ 1), 
the cooling is negligible, and the compressed parcel be- 
haves as a quasi-adiabatic sound wave. This case is known 
to be unstable only if the adiabatic criterion is satisfied, 
which for our piecewise power-law cooling function re- 
quires /3 < 1/(1—7) = —3/2. In the opposite case {r] <C 1), 
the PdV heating is negligible in comparison with the cool- 
ing, and the pressure is determined by the thermal equi- 
librium condition. In this limit, the velocity perturbations 
are unstable when the isobaric condition is satisfied, be- 
cause they behave as a condensation mode (see, e.g., Shu 
1992; Paper III). Thus, velocity fluctuations strong enough 
that they dominate rj are linearly stable at ry > 1 (small 
scales) but unstable at ry < 1 (large scales) in the atomic 
ISM. Moreover, when Tu < Tg, the instability of large-scale 
perturbations manifests itself in a null resistance of the 
flow to turbulent compressions rather than in spontaneous 
growth of the perturbations, because the bulk velocity is 
larger than the sound speed, and so the driver of the den- 
sity growth is the turbulent velocity compression rather 
than the thermal pressure gradient. 

It is worth remarking that velocity fluctuations arc a 
natural way of producing the density fluctuations that 
constitute the initial conditions for the subsequent devel- 
opment of TI, because the gas obeys to the continuity 
equation (mass conservation), and thus the production of 
a density enhancement requires a converging velocity field. 

The case of strongly nonlinear compressions probably 
depends on the closeness of 77 to unity. Ifr]':^ 1, then even 



very strong compressions behave adiabatically, and arc not 
able to induce very large density enhancements. If rj is still 
larger than unity at the higher density, triggering the iso- 
baric mode will be extremely difficult. On the other hand, 
if T] is not much larger than unity initially, then the density 
increase induced by the compression may raise the cooling 
rate to large enough values that rj becomes smaller than 
unity, and triggering of the isobaric instability ensues. 

4. NUMERICAL RESULTS 

4.1. Evolution of a single density perturbation and final 

stages of TI 

In this section we attempt to quantify the dynamical 
nature of the condensation process in the context of the 
formation of interstellar clouds by TI. In particular, we 
seek to determine the physical perturbation sizes at which 
the condensation develops supersonic velocities, the times 
needed for these velocities to subside, and the state of the 
leftover diffuse gas, which, as we shall see, contains signif- 
icant fractions of mass still in the unstable density range. 

Let us consider first how the nonlinear evolution of a 
Gaussian density perturbation depends on r]o- To this end, 
we numerically follow the evolution of perturbed density 
fields of the form 

po{x) = l+Ao exp . (7) 

where the density is in units of cm~^. We take Aq = 0.2 
and consider three values of a such that the FWHM of 
the Gaussian profile is 3, 15 and 75 pc. We refer to the 
corresponding simulations as DEN3, DEN15 and DEN75. 
According to our discussion in §3.2, it is expected that 
for model DEN3, the evolution will be closer to isobaric 
and quasi-static (?7o ~ 3.7, asstuning Ts = FWHM/2cs), 
whereas dynamical compression regimes are expected for 
DEN15 (r/o « 0.75) and DEN75 {rjo w 0.15). AU these 
simulations start with zero velocity. 

Model DENS starts with constant pressure. The ini- 
tial pressure is P/k = 2500 K cm~^, which equals the 
equilibrium pressure Peq at 1 cm""^. Figure 2 shows the 
density, velocity and pressure profiles at four different tem- 
poral snapshots. The maximum of the velocity, fmax(i), 
increases continuously until t ~ 3.8 Myr, at which time 
^'max ~ 0.86 km/s (Mach number 0.18). As expected, al- 
though the pressure is not perfectly constant (see Fig. 2b), 
the evolution can be considered as quasi-isobaric, as can 
be seen in Figure 4. Indeed, the maximum pressure vari- 
ation is \6P\ /Po « 0.2. As a consequence, there is no 
appreciable density overshoot beyond pisob (§3.2), and the 
final pressure inside the cloud recovers a value very close 
to the ambient pressure. 

By time f » 4 Myr, a nearly stationary "slow-growth" 
stage is reached, in which most of the mass is in the cen- 
tral condensation at nearly zero velocity, with low density 
material falling in from the surrounding medium, but de- 
celerating smoothly as it enters the cloud, without the for- 
mation of shock. The time for reaching this stage increases 
significantly when the initial perturbation amplitude is de- 
creased (^ 5.5 Myr for Aq = 0.1). In order to quantify the 
amount of mass within the condensation, we integrate the 
density over the length of the region where p > lOpo, lead- 
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ing to a column density^ N ~ 0.045 M© pc~^ at t = 4.2 
Myr. This value is 67% of the column density initially 
within \x\ < a. 

It is important to remark that in the "final" slow-growth 
state mentioned above, the leftover diffuse gas is far from 
static, with velocities of up to 0.7 km s~^. Although the 
accretion rate onto the condensation is low because of the 
low density of the accreting gas, the amount of mass con- 
tained in this regime is still large, ~ 8 times larger than the 
amount of mass in the condensation. This result is in part 
an artifact of the fact that we considered a single pertur- 
bation, which by itself does not consume a large amount of 
mass by the time it reaches pisob- Smaller final fractions of 
mass still in the unstable range are expected in the case of 
multiple perturbations, but even in this case, much longer 
times are needed to evacuate the unstable range (§4.2). 
Smaller unstable fractions are also expected if the back- 
ground medium is already in the stable diffuse phase and 
the condensation is triggered by a strong compression (we 
thank P. Hennebelle for pointing this out). 

Another point worth remarking is that the accreting gas 
evolves in a nearly isobaric regime in response to the insta- 
bility, but such that the small existing pressure variations 
have the same sign as the density variations (Fig. 2), and 
thus this gas has no further tendency to fragment, even 
though its density lies in the "unstable" range. Moreover, 
this implies that this density range does not precisely cor- 
respond to the unstable temperature range, as shown in 
Fig. 3. 

Let us now consider the case of models DEN15 and 
DEN75. In these, we set up the initial pressure to the 
corresponding foq values, Pq = P^q {po{x)), which is rea- 
sonable for models with low /y's, i.e. short cooling times. 
Figure 5 displays four snapshots of DEN15. We can distin- 
guish three evolutionary stages: (1) An initial compression 
stage in which all the parcels in the flow are in the ther- 
mally unstable range, the fluid is continuously accelerated, 
and the pressure gradient becomes increasingly steeper, 
eventually generating transonic infall motions. (2) A re- 
pressurization, or "crushing" , stage, starting when the core 
of the cloud enters the thermally stable temperature range 
and the pressure starts to increase in the center of the 
cloud, while the outside material is still unstable and con- 
tinues to accelerate inwards, generating a shock at the edge 
of the cloud at the time of maximum compression. Note 
that this is a radiative shock occurring over distances of 
the order of the cooling length, and so it is well resolved by 
the numerical grid. (3) A slow, subsonic accretion stage, 
which slowly fades away. During the initial compression 
stage the system is closer to thermal equilibrium than in 
DEN3 (see Fig. 6). The maximum Mach number achieved 
in this run was 0.6, and occurred at 6.1 Myr. It is neces- 
sary to start with a perturbation of FWHM of 40 pc for 
the infall velocity to be transonic. 

The shock wave that emanates from the cloud bound- 
ary is worth discussing. This shock is rather weak, as can 
be seen from the pressure jump in Figure 5b, but has im- 
portant effects in the warm medium. As shown in the 
P-p diagram (Fig. 6), as the shock propagates outwards, 
it heats the inflowing gas, bringing it closer to an isobaric 



regime, similar to the one of run DENS, while bringing it 
out of thermal equilibrium. In fact, the relative pressure 
variation SP/Pq from x = to the position just behind the 
shock is only ~ 0.2 in the slow accretion stage, whereas the 
temperature varies from 50 to 3000 K over this same in- 
terval. The column density of the cloud is 0.26 Mq pc~^, 
about 76% of the initial column density within \x\ < a 
(0.335 M0 pc~^). Since the ram pressure in the bound- 
ary of the cloud is larger than in model DEN3, the cloud 
density is slightly larger too. 

We now consider run DEN75, for which rjo ~ 0.15. Fig- 
ure 7 shows the evolution of density, pressure and velocity. 
The three evolutionary stages discussed in model DEN15 
are also present in DEN75. In this run, the gas behaves 
precisely along the thermal equilibrium P-p curve in the 
initial collapsing stage {t <21 Myr), as we can see in Fig- 
ure 8. In this situation, the gas acquires the maximum pos- 
sible values of the infall velocity imposed by the develop- 
ment of pure TI. In fact, at the time when the cloud begins 
to reach the thermally stable phase (T « 140 K), the max- 
imum velocity is 3.0 km s~^, and the maximum local Mach 
number is 1.2. The maximum of the velocity continues in- 
creasing up to 3.6 km s~^, agreeing with our estimate for 
■^inst in §3.2 to within 10%. Consequently, the supersonic 
crushing process is more violent than in DEN15, causing 
a large density overshoot, of roughly 30 x pisob- Note that 
the actual value of the density peak in these simulations 
is limited by slioc;k viscosity and mass diffusion, although 
this value is of little interest in ID contractions because, 
as mentioned in §3.2, it is not viable for triggering gravi- 
tational instabilities, except, perhaps, if the densities are 
high enough to trigger a non-equilibrium instability in the 
cold gas due to molecular cooling (Koyama & Inutsuka 
2000), which we do not include here. The density of the 
cloud will subsequently decrease slowly with time, but it 
will remain significantly larger than pisob (sec below). 

During the compression phase, the small-scale Fourier 
components of the Gaussian perturbation blow up due to 
their larger growth rates, forming two (symmetric) promi- 
nent density peaks, which ultimately merge with each 
other (Fig. 7; note that only one peak is seen, since only 
one side of the simulation is shown) . Their formation de- 
pends on details of the initial conditions, being nonexistent 
in the presence of a small initial convergent velocity field. 
The overall evolution of a simulation (DEN75b hereafter) 
with FWIIM=75 pc, Aq = 0.7 and initial velocity field 
vq oc —VPq/po, with initial maximum Mach number of 
0.4, does not present remarkable differences except for the 
presence of stronger accretion. 

As with run DEN15, at later times, the shock wave prop- 
agates outwards from the cloud boundary, causing a tem- 
perature jump of AT ~ 300 K across it, nearly restoring 
pressure equilibrium, albeit at the expense of throwing the 
outside gas out of thermal equilibrium. Figure 9a shows 
the density, the bulk and sound speeds, and the cooling 
times out to 75 pc from the center at 27.0 Myr. Within 
the region 1 < j.xj < 75 pc from the center, the density lies 
in the range 0.5 < p < 3 cm~^, with cooling times between 
0.8 and 2.7 Myr. Within 20 < < 75 pc, the cooHng time 
is roughly constant (tcooI ~ 2.0-2.7 Myr), and the typical 



^ Note that the units of column density are Mq pc ^ even though the simulations are one-dimensional, because the physical units of the 
density are g cm~^ anyway. 
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value of T] with respect to the size of this region (55 pc) 
is ~ 0.8-0.9. This value of 77 explains the approximately 
isobaric evolution of the flow across this region, which is 
characterized by relatively low densities (p < 0.7 cm~^). 
At intermediate densities, i.e. within 3 < |a;| < 20 pc, 77 
also has values ~ 1 (again with respect of the region size, 
~ 17 pc), because both TcooI and Tg decrease towards the 
cloud. Thus, this simulation has evolved towards values of 
Tj closer to unity over most of the medium external to the 
cloud, except at the cloud core (|a;| < 0.25 pc), where the 
condition 77 <C 1 is preserved (§3.2). Note also that for 
perturbations significantly smaller than the sizes of the 
regions mentioned above, r] is actually ^ 1 in the warm 
medium. 

The stage of slow accretion at the cloiid boundary is es- 
tablished after ~ 25 Myr. In Figure 9b the central part of 
the simulation, including the accretion front, is shown at 
t = 27.0 Myr. It is worthwhile to note that there is still 
an appreciable excess of density with respect to pisob at 
that time (p > 2pisob)- The column density of the cloud is 
1.65 Mq pc~^ at 27.0 Myr, very similar to the initial col- 
umn density within |a;| < a (1.68 Mq pc~^), whereas the 
mass column density remaining in unstable regions with 
p < 3 cm-3 and 1400 <T < 6000 K (not including the 
region unaffected by the perturbation) is 3.2 Mq pc~^. 
In fact, only a small fraction of the low-density gas has 
reached densities corresponding to the warm stable phase 
(Fig. 7) (no gas at all reached that phase in runs DEN3 and 
DEN15). As mentioned with respect to run DEN3, this is 
in part a consequence of having considered a single density 
perturbation in an unstable, homogeneous medium. 

The final slow-growth state of this large-scale pertur- 
bation is actually very similar to that of DEN3, after the 
shock has established a quasi-isobaric regime. The gas still 
accreting has subsonic velocities, is out of thermal equilib- 
rium, and is traversing the "unstable" density range, al- 
though with density and pressure fluctuations of the same 
sign, so that it has no further tendency to fragment. Due 
to the low accretion rate, long times will be required to 
evacuate the unstable range, as already pointed out by 
Goldsmith (1970). 



4.1.1. 



Convergence study 

In order to assess the effect of numcirical resolution, we 
have conducted three realizations of each model DEN3, 
DEN15 and DEN75b, in which the resolutions differed by 
factors of two^. For DEN3, the number of grid points, Ux, 
were 800, 1600 and 3200, which correspond to a resolution 
around a; = of 3 x lO^^, 1.5 x 10"^ and 7.6 x lO"-"^ pc, 
respectively. Almost identical results were found (differ- 
ences of less than 1% in all variables). For DEN15 we used 
1800, 3600 and 7200 grid points, with associated resolu- 
tions 7.6 x 10-2, 3.8 X 10-2 and 1.9 x 10-^ pc, respec- 
tively. Again, the variations were insignificant (less than 
1%). For DEN75b, simulations were conducted with 2300, 
4600 and 9200 grid points. In Fig. 10 we compare two 
snapshots for = 2300 and rix = 9200 at an advanced 
time when steep gradients are already present. The den- 
sity and velocity profiles are similar in both cases, whereas 
the pressure is slightly less smooth in the low resolution 
simulation. In the high-resolution simulation the pressure 

* We present here the case of DEN75b because it is the most computationally demanding situation. 



is better resolved but the velocity field is unaffected. 

It is possible that the growth of small-size structures 
consisting of a few grid points are limited by resolution. 
To be sure that the results are not sensitive to resolution, 
experiments with increasing resolution have been carried 
out for all of the simulations presented throughout this 
paper. High enough resolution has been used to ensure 
that it does not affect the results. See also Paper III for 
further discussion. 

4.2. Multi-scale density perturbations 

As mentioned above, interesting phenomena may be ex- 
pected if the initial density profile contains a spectrum 
of wavelengths (see also Burkert & Lin 2000). Small scale 
perturbations (77 > 1 ) have larger growth rates than larger- 
scale ones (77 < 1) and thus will quickly outgrow them. 
Suppose now that we have two separate scales, one with 
T] and the other one » 1. If the clouds formed from 
small scales do not consume an important fraction of the 
mass, their bulk motions will be dominated by the pres- 
sure structure of the larger scales. As a consequence, these 
clouds will behave as if passively advected by the larger 
ones and will be accelerated towards pressure minima. In 
this case the evolution may be envisaged as a "linear" su- 
perposition of the evolution of small and large scales. In 
a more general situation, however, the evolution could be- 
come extremely complex because of the competition be- 
tween the growth of perturbations of different sizes. To 
quantify the relevance of this competition on the dynam- 
ics and mass spectrum of the forming clouds, we consider 
a superposition of fluctuations with wavenumbers in the 
range fcmin < k < /Cmax, where kmin = 27r/Z, with I the box 
size, and fcmax is varied to consider various extents of the 
perturbed range of scales. The amplitudes of the fluctua- 
tions follow a power-law of the form pk = Co(fc/fcmin)~*- 
We consider two values of the ratio fcmax/fcmin and two val- 
ues of the spectral slope S. Moreover, we consider another 
set of simulations in a box of size 10 times smaller. Note 
that, if the amplitude of the initial perturbations is an in- 
creasing function of the wavenumber (i.e. S < 0), it is clear 
that the small scales will be able to consume an important 
fraction of the gas with no chance for large scale perturba- 
tions to survive. Thus, we restrict our discussion to (5 > 0, 
in particular to the values 3/2 and 1/2. For these values 
of 6 we now discuss the conditions under which the small- 
scale perturbations do not appreciably affect the overall 
evolution. 

First, we compare the evolution of two runs with S ~ 
3/2, labeled SP20a and SP4a (for "superposition"; sec Ta- 
ble 1). These runs have ratios kmax/kmin = 20 and 4, 
respectively. We consider a box size of 250 pc, which im- 
plies that, for these runs, 7;o(A:min) — 0.04. Four snapshots 
of their evolution are displayed in Figure 11. The num- 
ber of condensations at intermediate times is 16 and 5 for 
SP20a and SP4a. respectively. This number depends al- 
most linearly on fcmax/fcmin for these simulations. 

The problem of competition for the available mass can 
be summarized in the question of whether clouds formed 
from small scales in run SP20a group together to form a 
density pattern that at large scales is similar to that of 
run SP4a. If that happens, we can say that perturbations 
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at small scales arc unable to suppress the development 
of TI at larger scales. Note that we speak in terms of 
Fourier wavelengths only for the initial spectrum. In con- 
trast, once condensations are growing, we focus on their 
clustering and distribution in physical space. 

From panels (c) and (e) of Figure 11, we see that some 
of the incipient condensations in run SP20a are disrupted 
when they are embedded in a diverging flow (for instance 
the clumps between x = —25 and x = 20). Conversely, 
structures placed in a compression zone grow very fast. 
Note that some of the clouds in the group near x = 60 
are moving relatively fast (~ 1.5 km/s) at 25 Myr. At 
later times, the clouds undergo a series of mergers that 
ultimately lead to four main condensations at positions 
X = -110, X = -95, -40 and 60 at f « 38 Myr. 

For this run, the column density ratios of gas mass 
between the warm, the thermally unstable and the cold 
phases are 1 : 0.5 : 1 at 15 Myr and quickly evolve to 
0.69 : 0.025 : 1 at 20 Myr, whereas at 25 Myr they are 
0.46 : 0.01 : 1. 

As expected from their larger sizes, the growth rate of 
the density structures is smaller in run SP4a, but each 
condensation accretes more matter. In fact, most of the 
gas is still in the thermally unstable phase at 15 Myr, with 
column density ratios 0.9 : 9.7 : 1. At times ^ 25 Myr, 
these ratios are 0.66 : 0.034 : 1. Still, the fraction of gas in 
the thermally unstable phase is higher than in SP20 at 20 
Myr. At those times 25 Myr), four main condensations 
have already formed at a; = —100, —40, 35 and 70. Com- 
paring with SP4a, the evolution of run SP20a is somewhat 
more dynamical, in the sense that clouds acquire larger ve- 
locities and are prone to suffer mergers more easily than 
in run SP4a, in which clouds are more massive and static. 
The locations of the condensations are the same except 
for those at a; = 35 and x = 70 which have not had time 
to merge yet; they will do so in an exceedingly long time, 
at < w 90 Myr. We conclude that the final states are not 
equivalent in this case. 

Next, we consider two more sets of simulations using 
the same density spectrum, but varying the simulation 
box size. When the box is five times larger, the differ- 
ences are even more significant. However, if the box size 
is reduced a factor ten (25 pc, runs SP20b and SP4b, re- 
spectively), so that for ry(A:min) = 0.4, the final states of 
the two initial spectra are much more similar, as seen es- 
pecially from their velocity fields (see Fig. 12). Hence, the 
evolution is not only sensitive to the initial spectrum of 
the perturbations, but also to r?(A;min)- 

Let us now consider a set of similar experiments but with 
^ = 1/2. First, we take ?7o(fcmin) = 0.04 (runs SP20c and 
SP4c). In Figure 13 we see that nm SP20c contains sev- 
eral condensations that cannot be associated to any main 
cloud in SP4c. The most noticeable difference between 
SP20c and SP4c is that SP20c is able to group condensa- 
tions into two clusters whereas the four condensations in 
SP4c would need an exceedingly long time to coalesce, as 
can be seen by comparing the velocity gradients between 
positions x = —110 pc and x = —50 pc for the two runs 
(Figure 13). Hence, the final states of runs SP20c and 
SP4c are qualitatively very different. 

Next, we consider a set of runs with 6 = 1/2 and 
??o(fcmin) = 0.4 (runs SP20d and SP4d) and another set 



with r7o(fcmin) 4.0 (runs as SP20e and SP4e). Run 
SP20d still exhibits noticeable differences with respect to 
SP4d. In Figure 13 we see that there are three condensa- 
tions between x = 1 and x = 7 pc, generated from small 
scale perturbations in SP20d, but they are absent in SP4d. 
These condensations have positive velocities, whilst the 
gas between x = 1 and x = 5 has the opposite veloc- 
ity in SP4d. This implies that in this case the emergence 
of small-scale condensations produces important modifica- 
tions on the evolution of the gas at large scales because the 
late-time mass distribution is significantly altered if they 
are present. However, the opposite conclusion is reached 
for SP20e and SP4e, which are seen to have more simi- 
lar states, again judging in particular from their velocity 
fields (see Fig. 13). The density peaks between x — and 
X = 0.75 pc in run SP20e do not grow at all between 6 
and 8 Myr. 

We conclude that both the spectrum exponent and 
'7(fcmin) determine the outcome of a superposition of modes 
of various scales. It appears that, as the initial density 
power spectrum is shallower, the perturbations need to 
extend to smaller scales (larger 77's) in order for the small- 
scale fluctuations not to strongly disturb the growth of 
the larger-scale ones. In particular, we have found that, 
for (5 = 3/2, a value r]{kinin) = 0.4 is required, while for 
(5 = 1/2, a more stringent condition, ry(fci„in) = 4.0 is nec- 
essary. These critical values could be different for 2D or 
3D systcinis. 

It is worthwhile to point out that even though the devel- 
opment of TI may be suppressed at large scales due to mass 
consumption by small scales, massive clouds may form by 
coalescence of small clouds (see e.g., Elmegrecn 1990, Mur- 
ray & Lin 1996). As mentioned before, clouds formed at 
small scales may present appreciable bulk velocities that 
facilitate the merger of clouds. This effect has recently 
been seen also in 2- and 3D simulations by Koyama & In- 
utsuka (2002) and Kritsuk & Norman (2002). The coales- 
cence is exaggerated in one dimensional systems, because 
the merger rate does not depend on the filling factor of the 
clouds. However, in more than ID, the merger rate may 
still be large if the small-scale clouds are advected towards 
the condensation center of the large-scale ones (see §5.2). 
In general, our ID simulations show that small scale per- 
turbations do not lead to a "forest" of small clouds because 
they coalesce, building more massive ones. 

4.3. Simulations of forced flows 

We now consider the case of flows forced continuously 
in time, representative of a turbulent ISM, to determine 
whether turbulence may oppose the condensation process 
(cf., §3.3) under some circumstances, and, if so, at which 
scales and values of the rms Mach number. In a domain 
of size I, we force the flow at scales large compared to 
the final sizes of the condensations. This implies that, 
once formed, the condensations cannot be disrupted by 
the random forcing. Note, however, that in higher dimen- 
sionalities dynamical instabilities could disrupt the clouds 
while they form, or after they have formed; see §5.2 for 
further discussion. 

We adopt a random force / of the form 

f{x,t) = m{Nexp[ik{t)x + i(j){t)]}, (8) 
where 3? denotes the real part, k{t) is a time dependent 
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wavenumber, (j)(t) is the phase and x the position. Follow- 

1/2 

ing Brandenburg (2001), we take N = foCg [k{t)cs/6t] ' , 
where /o is a constant factor and St is the length of the 
timestep. The values of A;/(27r//), where / is the box length, 
and of (p are selected at each timestep randomly in the 
ranges [3,10] and [0,27r], respectively. This implies that 
the scale of maximum forcing is At = I/IO. The positive 
exponent (1/2) in TV implies that the maximum strength 
of the forcing occurs at the highest wavenumbers fc, i.e., 
roughly at 1/10 of the box size. A summary of the forced 
simulations is given in Table 2. There, Tcond is defined as 
the time for a condensation to reach a density p = lOpo 
for the first time, (e) is the ratio of the time-averaged in- 
ternal energy (with the average performed over Tcond) to 
its initial value. A^rms is the rms Mach number, also av- 
eraged over Tcond to increase statistical confidence. Note 
that Alrms fluctuates strongly in time, in contrast to e, 
which behaves much more smoothly. The different values 
of Tcond in one row in Table 2 correspond to runs with 
different seeds for the random numbers. 

Let us first consider experiments with I = 3 pc, in which 
Af = leq/30, since ^oq ~ 10 pc (cf. §3.2). For low Mach 
numbers (say Ai < 0.1), it then holds that the sound 
speed is larger than the typical values of the velocity fluc- 
tuations, with 77 ^ 1. Consequently, the condensations be- 
have essentially as entropy perturbations, and evolve ap- 
proximately isobarically (cf. §3). Interestingly, however, 
the various realizations of run 1, which have I = 3 pc 
and M. = 0.03, present a mean Tcond ~ 12 Myr, about 
three times the condensation time for the case with an 
initial density perturbation of amplitude 20 % and the 
same size. A large dispersion of Tcond is apparent for flows 
with A^rms ~ 0.15-0.3, suggesting that the formation of 
condensations is sensitive to the past history of the flow, 
becoming somewhat fortuitous for higher Mach numbers. 
For even larger values of A^rms (^ 0.3), eventually the 
condensation is suppressed altogether (runs 5-9). In these 
cases, the internal energy initially increases with time un- 
til it starts to fluctuate 4%) around a value more than 
twice its initial value. This occurs when cooling is able to 
compensate the thermal energy injection by the random 
forces, since the net cooling rate increases also with time 
as both the strength of the perturbations in density and 
temperature increase. Specifically, energy injection and 
net cooling are in balance with a mean value of 0.038 erg 
g~^ in run 7 after a time ^ 0.8 Myr, whereas this oc- 
curs after a time ^ 0.35 Myr and for a value of 0.175 erg 
s~^ g~^ for run 9. Our results thus suggest empirically 
that, for the cooling functions we use, the transition to 
stabilization occurs roughly at c ^ 1 (Airms ^ 

0.3). 

On the other hand, for the case I = 100 pc (run 17), it is 
clear that the forcing is applied at scales where ?? < 1, and 
indeed the density fluctuations generated by the random 
compressions rapidly cool and acquire a temperature and 
a pressure close the thermal equilibrium values. Since the 
average density is in the unstable range, TI proceeds unim- 
peded, manifesting itself not so much as a tendency for 
the fluctuations to grow, but in the null resistance of the 
flow to external compression (§3.3). Indeed, the random 
compressions promote the instability, as indicated by the 
decrease in Tcond with increasing AJrms- This is the regime 
of the randomly forced simulations of Paper I, and corre- 



sponds to the case of an effective imaginary sound speed 
of the wave mode, which consequently does not propagate 
(e.g., Field 1965; Shu 1992). 

The case with Z = 10 pc (runs 10-16) has Af ^ Zcq/10. 
Nevertheless, Table 2 shows that the condensation of the 
density fluctuations cannot be inhibited anyway, even for 
M 

rms 0.6. We see that e can hardly be increased by 
factors of 1.9 in this case. We understand this as a conse- 
quence of the fact that, for such large Mach numbers, the 
reduction of the cooling time by the compression-driven 
density enhancement can reduce rj, allowing the triggering 
of the instability under the isobaric mode by the nonlinear 
velocity fluctuations. It appears that this does not happen 
in the case ^ = 3 at large Alrms because rj is large enough 
in that case that the shocks produced are too close to adi- 
abatic and cannot raise the density enough to destroy the 
adiabaticity. We should remark that for simulations with 
^ = 10 pc and strong forcing (Afims > 0.5, e.g. nm 16) 
compressions can produce transient density fluctuations 
~ lOpO) which are dispersed by their own pressure. 

We conclude that, as conjectured in §3.3, the growth of 
perturbations can be suppressed provided that the forcing 
is applied at scales small enough that 77 ^ 1 and produces 
a large enough rms Mach number so that r] is dominated 
by the turbulent velocity fluctuations rather than by the 
sound speed. This is because, under these conditions, the 
perturbations are nearly adiabatic, and are (at least) lin- 
early stable. The fact that even large rms Mach numbers 
in runs with ^ = 3 pc do not induce condensations indi- 
cates that the stability of these perturbations goes well 
beyond the first order. 

5. DISCUSSION 
5.1. Comparison with previous work 

Some of the issues considered in this paper have been 
discussed by previous authors in the context of the frag- 
mentation of proto-galactic and proto-globular cluster 
clouds (e.g., Murray & Lin 1989, 1990), cooling fiows in 
galaxy clusters (e.g., David et al. 1988; Brinkmann et al. 
1990; MalagoH et al. 1990; Kang et al. 1990) and the ISM 
(e.g., Hennebelle & Perault 1999, 2000; Burkert & Lin 
2000; Koyama & Inutsuka 2000). In this section we put 
our results in perspective with respect to those works. In 
Paper I we have presented a comparison between our ISM- 
oriented scenario and previous work considering quasi- 
static development of TI. 

It should be first pointed out that some of those stud- 
ies have solved the chemical rate equations, allowing for 
out-of-chemical-equilibrium abundances of the relevant 
coolants (e.g., Murray & Lin 1989, 1990; David et al. 
1988; Kang et al. 1990; Koyama & Inutsuka 2000). This is 
crucial in situations in which considering out-of-chemical- 
equilibrium cooling rates can make the difference between 
stability and instability. However, it is not so crucial in 
relation to the present-day ISM, which is already generally 
agreed to be thermally unstable in the range 300 K < T < 
6000 K under chemical equilibrium conditions. Once the 
medium is known to be unstable, it seems reasonable to 
think that variations in the cooling rates due to devia- 
tions from chemical equilibrium should only cause mod- 
ifications in the details of the condensation process, but 
that the main qualitative behavior should not change sig- 
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nificantly. For example, the development of supersonic 
motions by the development of TI similar to our results in 
§4.1 in the context of cooling flows and of proto-globular 
cluster clouds has been observed by David et al. (1988) 
and Kang et al. (2000), respectively. Moreover, Murray & 
Lin (1989, 1990) discussed the effects of various types of 
perturbations, cases with and without background heat- 
ing, and the need for a minimum cloud collision veloc- 
ity (equivalent to a compression) in order to trigger TI 
in proto-globular cluster clouds subject to external heat- 
ing. The latter result is similar to that of the more recent 
study by Hennebelle & Perault (1999) in the context of the 
ISM, even though these authors did not explicitly solve for 
the chemical rates. Thus, similar phenomenology is found 
with or without non-equilibrium chemistry, once the sys- 
tem is unstable. In any case, the caveat of our neglect 
of non-equilibrium chemistry and rates should be kept in 
mind. 

Now let us specifically consider recent studies of TI in 
the context of the ISM. First, as mentioned above, the 
studies of Hennebelle & Perault (1999) and Koyama & 
Inutsuka (2000) were again concerned with the trigger- 
ing of TI in the stable phases of the ISM by means of 
strong shocks or compressions. These studies are related 
to the issues we discuss in §§3.3 and 4.3, although in our 
case the medium is already isobarically unstable from the 
start. Interestingly, however, the case of velocity pertur- 
bations with f] > I is somewhat resemblant of the case 
of an initially stable medium, in the sense that a strong 
enough compression is needed in order to raise the cool- 
ing rate to high enough values that condensation can oc- 
cur. Also, it is possible that the induction of TI in the 
cold neutral medium due to molecular cooling found by 
Koyama & Inutsuka could occur during the transient over- 
shoot we found in the condensation of large-scale perturba- 
tions if we included such cooling, but in any case this phe- 
nomenon would only support our suggestion that the static 
pressure-equilibrium configuration is hard to achieve, since 
the diffuse clouds could be subject to a further instability 
that would transform them all the way to small molecular 
cloudlets. 

On the other hand, Burkert & Lin (2000) did aim their 
study at the formation of clumpy clouds in an unstable 
ISM, but they considered a different physical situation 
than the one we have considered here, a fact which caused 
them to reach significantly different conclusions from us. 
Specifically, they considered a case with no background 
heating, a choice that changes the instability criterion and, 
more importantly, implies that the whole flow is cooling, 
so that the development of TI consists of a perturbation 
cooling faster than the background medium. Thus, per- 
turbations that cool isochorically and grow slowly (i.e., 
large-scale ones) cannot get "too far ahead" of the back- 
ground gas before a minimum temperature is reached at 
which the gas either ceases to cool or exits the unstable 
range. As a consequence, Burkert & Lin concluded, simi- 
larly to Murray & Lin (1990) in the context of protogalac- 
tic clouds, that such perturbations cannot reach very large 
density enhancements over the background, and that only 
isobarically cooling ones do. This is in sharp contrast to 
our conclusion that it is precisely the large-scale fluctua- 
tions that reach the highest amplitudes (albeit only tran- 



siently). The reason is that in our case the background 
is not cooling, and thus there is no competition between 
it and the fluctuations to reach the minimum tempera- 
ture. Which case applies at a given location in the real 
ISM should depend on the presence or absence of a pop- 
ulation of heating sources that maintain a relatively con- 
stant background heating rate. It is in general agreed that 
the ISM is permeated by a roughly uniform UV field and 
cosmic rays that provide a background heating, although 
non-negligible variations are certainly expected to occur 
(Hollenbach, Parravano & McKee 2001). 

5.2. Adequacy of the one-dimensional approach and 

two-dimensional simulations 

The one-dimensional numerical approach that we have 
used here clearly has some limitations. It was chosen 
in order to fully resolve the condensations even at the 
time of maximum compression, but it has the disad- 
vantage of being unable to capture effects that involve 
solenoidal (vortical) motions in the flow. In particular, 
it has been suggested that dynamical instabilities such 
as Raylcigh- Taylor (RT), Ritchmyer-Meshkov (RM), or 
Kelvin-Hclmholtz (KH) may destroy the clouds during or 
after their formation in stratified media (e.g., Malagoli et 
al. 1990; Hattori & Habe 1990; Reale et al. 1991; Murray 
et al. 1993; but see Vietri et al. 1997 for an opposite view). 
However, within the scope of the present paper, in which 
we have chosen to study the development of TI in isolation, 
we are not concerned about the future fate of the clouds 
formed, but only about issues of the cloud formation pro- 
cess within the idealized scenario that TI is responsible for 
this process. This amounts to studying the limiting, most 
favorable case in which the classic scenario of an equi- 
librium ISM, with thermal-pressure bounded clouds and 
sharp phase transition could arise, and our main conclu- 
sion is that even within this scenario such a picture of the 
ISM appears unlikely. 

Perhaps the closest point of contact of the issues dis- 
cussed in this paper with the possibility of cloud destruc- 
tion by dynamical instabilities occurs in the context of our 
discussion of the evolution of velocity perturbations. As 
we have seen, these may avoid condensation for arbitrarily 
long times if they have moderately large amplitudes and 
occur at small enough scales. Thus, in a sense, the velocity 
fluctuations even in our simple ID case have the effect of 
counteracting the development of the isobaric mode of TI 
analogous to that of the dynamical instabilities discussed 
by other authors. 

It is thus of interest to compare the effects of those dy- 
namical instabilities on the condensed clouds with the ef- 
fects of the random forcing we have considered here. Con- 
cerning the KH instability, it is expected to occur as the 
condensations begin to sink in a stratified medium as they 
develop inverse buoyancy (Balbus & Soker 1989). Thus, it 
could only start to act after the condensation has reached 
a sufficiently advanced evolutionary stage. That is, its ef- 
fect is one of destroying a cloud after it has formed, rather 
than obstructing its growth, as we have found here. 

On the other hand, the RT instability might indeed oc- 
cur while the condensation proceeds, as the flow is acceler- 
ated inwards of the condensation, with the density gradi- 
ent pointing also inwards. To test this possibility, we have 
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performed two two-dimensional simulations analogous to 
run DEN75, but including weak random initial velocity 
perturbations with rms Mach numbers 0.1 and 0.01 with 
a half- wavelength of 1/8 of the box. We find no traces of 
development of RT. Instead, the velocity fluctuations gen- 
erate density fluctuations that evolve in a very similar way 
to the simulations of §4.2, with the small-scale fluctuations 
condensing first, and then merging towards the center of 
the large-scale perturbation. The certainty that RT does 
not arise comes from the fact that the velocity field con- 
verges onto the fragments rather than shear between them, 
as would be expected for RT. Thus, we conjecture that TI 
development inhibits that of RT. 

Finally, the RM instability probably docs not apply to 
single clouds formed by TI because it occurs upon the in- 
teraction of a shock with a discontinuity. Instead, what 
we have found in this paper is that the evolution of a sin- 
gle perturbation generates a shock wave that propagates 
outwards from the condensation. It may be relevant, how- 
ever, in the case of multiple fluctuations, as in §4.2, in 
which the shock propagating out of one condensation may 
interact with another. However, this again refers to the 
future evolution of the condensations formed by TI, which 
is out of the scope of this paper. 

A final note is that, in any case, the effects of all these 
instabilities go in the same direction as the issues discussed 
in this paper: as they tend to destroy the condensations, 
they tend to restore a more continuum-like ISM, rather 
than one separated in discrete phases. 

5.3. Implications 

The resiilts of this paper suggest a number of implica- 
tions. All of them require further investigation, but the 
possibility of their realization in the ISM is of great inter- 
est: 

1. The outcome of the development of TI, even if it 
were the sole cloud-forming agent in the atomic ISM, is not 
necessarily the production of quiescent thermal pressure- 
bounded clouds in contact with the warm stable phase. 
Large-scale density perturbations require long times to 
condense (over 7 Myr for perturbations of size > 10 pc), 
during which they are in a highly dynamic state, far from 
pressure equilibrium. At later stages, pressure equilibrium 
is approximately restored, but at the expense of disrupting 
thermal equilibrium of the surrounding medium throiigh 
the shock propagating away from the condensation. This 
surrounding gas continues to accrete onto the condensa- 
tion for much longer times, providing a ram presure that 
keeps the cloud density above the pressure-equilibrium one 
by factors of up to 2. Moreover, the accreting gas is not 
in thermal equilibrium, but in a nearly isobaric regime, in 
which the (small) pressure fluctuations have the same sign 
as those of the density, and therefore the medium has no 
strong tendency to fragment any further. In particular, 
this non-equilibrium situation implies that the unstable 
density range as shown in Fig. 1 does not exactly coincide 
with the unstable temperature range. 

Smaller-scale (large-vy) fluctuations (< 3 pc at (p) ^ 1 
cm""^, or even larger at proportionally lower densities) 
can condense and reach near thermal-pressure balance in 
shorter times (< 4 Myr) if they are nearly isobaric in na- 
ture, but they are stable if they are nearly adiabatic, a 



condition that can be accomplished if they originate from 
velocity fluctuations with sufficiently large A4rms at suffi- 
ciently small scales (§4.3). Moreover, small-scale fluctua- 
tions, if part of a full spectrum of fluctuations, may con- 
tinue to form part of the condensation of large-scale ones, 
if the relative amplitude of the former is small enough and 
they occur at sufiiciently large rj (§4.2). In general, it ap- 
pears that the pure nonlinear development of TI, even in 
the absence of other agents like magnetic flelds, rotation, 
stellar energy injection, etc., docs not necessarily lead to 
rapid formation of clouds with sharp phase transition at 
their boundaries. 

2. The existence of accreting gas in the "unstable" range 
at late stages of the evolution, and the stabilization pro- 
duced by velocity fluctuations in a certain range of scales 
suggests that in the turbulent ISM the "unstable" phase 
may actually be signiflcantly populated, since the density 
fluctuations are in general expected to have a dynami- 
cal origin. This may explain observational (e.g., Dickey, 
Salpeter & Terzian 1977; Kalberla et al. 1985, Spitzer & 
Fitzpatrick 1995, Fitzpatrick & Spitzer 1997; Heiles 2001) 
and recent numerical results by Gazol et al. (2001) and 
Kritsuk & Norman (2002), suggesting that relatively large 
amounts of gas are present in the ISM in the unstable 
temperature range. 

3. In summary, both the gas accreting onto the conden- 
sations at late times and the "unstable" gas in the presence 
of stable adiabatic fluctuations, are out of thermal equilib- 
rium. In general, it can be said that thermal and pressure 
equilibria are mutually exclusive for gas in the unstable 
range, because each implies a different value of the pres- 
sure at densities intermediate between po and pisob (fig- 
3). This incompatibility between the two equilibria can 
only be avoided by means of sharp phase segregation, i.e., 
when no gas in the unstable exists. However, if such gas 
exists, it should in general be expected to be out of ther- 
mal equilibrium. This suggests a possible observational 
test for determining whether the gas apparently seen at 
unstable temperatures corresponds to this kind of regime, 
or at least whether it is in thermal equilibrum or not. This 
could be done by either a) simultaneously detcrming two 
of its thermodynamic variables, or b) comparing directly 
observed cooling rates (e.g., fine structure lines) with the- 
oretical estimates of the heating rate (e.g., photoelectric 
heating) in specific regions (C. Heiles, J. Scalo, private 
communications). If this were confirmed, it would pro- 
vide strong evidence that turbulent motions populate all 
regions of the thermodynamic variable space, preventing a 
sharp segregation of the atomic ISM into the stable phases 
of TI. 

4. The nearly adiabatic response of the gas to velocity 
perturbations in cases with relatively large values of t], sug- 
gests that the warm diffuse ISM, in which this condition is 
satisfied with respect to perturbations of up to several pc 
(§4.1), may exhibit a very weakly compressible behavior, 
and thus be close to a Kolmogorov regime, which applies to 
incompressible turbulent flows. This may explain why ob- 
servations of the diffuse gas tend to find density fluctuation 
spectra with the signature of a Kolmogorov regime (e.g., 
Minter & Spangler 1996; Stanimirovic & Lazarian 2001; 
Dickey et al. 2001). The numerical results by Kritsuk & 
Norman (2002) also suggest a nearly adiabatic regime for 
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the diffuse gas. 

6. SUMMARY 

In this paper we have performed a detailed numerical 
investigation of the nonlinear development of the isobaric 
mode of TI in the context of the atomic ISM, and of its 
interplay with a spectrum of fluctuations in density or ve- 
locity, aimed at establishing the feasibility of the classical 
scenario of a thermal- and pressure-equilibrium ISM, in 
which the clouds are directly confined by the more dilute, 
stable warm gas, as suggested in the classical two- and 
three-phase models of the ISM (Field et al. 1969; McKee 
& Ostriker 1977), even under the idealization that TI and 
a spectrum of density or velocity fluctuations are the only 
processes at play. Most of the discussion is based on one 
parameter, rj, the ratio of the cooling to the dynamical 
times. 

We have found that the evolution for initial Gaussian 
perturbations of FWHM smaller than ^ 15 pc is approx- 
imately isobaric. In that case, a slow-growth stage char- 
acterized by quasi-stationary accretion of low-density gas 
is reached in « 4-6 Myr for initial amplitudes of 10- 
20 % over the mean density, although infall velocities 
> 0.7 km s^^ still exist at those times. The infall veloc- 
ities remain subsonic throughout the evolution for small 
perturbation sizes, but infall velocities with Mach num- 
ber > 0.5 are acquired if initial perturbations are large 
enough (> 15 pc). For the latter type of perturbations, a 
shock wave, formed during the "crushing" and overshoot- 
ing stage, propagates outwards from the cloud and throws 
the surrounding medium out of thermal equilibrium, es- 
tablishing near pressure balance. For the case with a Gaus- 
sian density perturbation with FWHM= 75 pc. the final 
mass density of the cloud is more than twice pisob 

due to 

the ram pressure of the infalling gas. The time to reach 

the slow-accretion stage ranges from 4 to 30 Myr for per- 
turbations with initial FWHM of 3 to 75 pc. These times 
are long enough that a major disturbance is likely to pass 
through the condensation. 

We have quantified the competition between scales by 
starting with a power-law spectrum of initial density per- 
turbations, to assess the possibility that the smaller scales, 
which have the fastest growth rates, could consume most 
of the mass available before large-scale perturbations have 
time to grow, thus establishing a "forest" of thermal- 
pressure bounded small clouds in relatively short times. 
We have found that small-scale perturbations may or may 
not use up most of the mass available depending on both 
the spectrum slope —5 and the initial values of rj for the 
smallest-scale perturbations. We have found that the in- 
clusion of large wavenumbers has a progressively smaller 
effect on the final state as 6 increases and as the pertur- 
bations are shifted to smaller scales, increasing r]. For 
example, for S = 3/2, a value of ri{k^in) > 0.4 is needed 
in order for the final states to be roughly equivalent, while 
for 6 = 1/2 this requires r/(/;',„i„) > 4.0. We looked at this 
problem to determine whether the dynamical character of 
the condensation of large-scale perturbations might be pre- 
served even in the presence of small-scale ones. Although 



our simulations indicate that the evolution is also very dy- 
namical even if small-scale perturbations are added, more 
work is needed to be conclusive concerning this result. 

However, we have found that a much more interesting 
effect appears in the case of a medium continuously stirred 
by an external force. We have argued that velocity fiuc- 
tuations are necessary to produce realistic density fluctu- 
ations in a continuum, their evolution thus ranging from a 
nearly adiabatic to a thermal equilibrium response as rj is 
varied from ^ 1 to ^ 1 . provided that in this case r/ is de- 
termined by the turbulent, rather than the sound crossing 
time. When the perturbations evolve nearly adiabatically, 
they cannot condense, because they re-expand before they 
can cool - i.e., they are stable. We have confirmed this 
suggestion numerically, showing that fiows with (p) = 1 
cm~^ forced at scales such that ry > 1 (A ~ 0.3-1 pc) and 
with A^rms > 0.3 do not develop any condensations over 
the entire duration of the simulations (between 12 and 26 
Myr). However, flows forced at scales such that < 1 
^ ^eq ~ 10 pc) produce condensations regardless of the 
forcing strength, because in this case either the effective 
waves have an imaginary sound speed (77 ^ 1), and trigger 
the condensation mode, which is unstable, or the density 
enhancement they induce causes rj to become smaller than 
unity by increasing the cooling rate (ry ^ 1). 

As implications, we have suggested that: 

1. The existence of accreting gas in the unstable den- 
sity range and the stabilization of the flow by small-scale 
velocity fluctuations may be at the origin of the relatively 
large amounts of gas mass in the unstable regime found in 
both observations (Dickey et al. 1978; Kalberla et al. 1985, 
Spitzer & Fitzpatrick 1995, Fitzpatrick & Spitzer 1997; 
Heiles 2001) and simulations (Gazol et al. 2001; Kritsuk 
k Norman 2002) of the ISM. 

2. Since transonic velocity fluctuations of scales up to 
several parsccs in the warm diffuse ISM have rj > 1, this 
medium should respond nearly adiabatically to them, and 
the flow regime should then be only weakly compressible. 
This may explain recent observational results suggesting 
that its kinetic energy spectrum is very close to the Kol- 
mogorov one (Minter & Spangler 1996; Dickey et al. 2001; 
Stanimirovic & Lazarian 2001), since the latter is expected 
to apply to incompressible or weakly compressible cases 
only. 

3. Simultaneous observational determinations of two 
thermodynamic variables should allow to distinguish 
whether the gas apparently seen in the unstable tempera- 
ture range is in thermal equilibrium or not, thus helping 
to decide whether equilibrium models of the ISM are ap- 
plicable or not. 
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Fig. 1. — Equilibrium curve for the adopted cooling function. We have labeled Pmax cLiid -Pmin the maxiniuiii a.iici miiiiniuni pressures 
the thermally unstable branch. The unperturbed density po, and the cloud density resulting from isobaxic evolution pisobi shown. 



16 




Fig. 2. — Time development of run DEN3. The density, pressure and velocity profiles at different times are plotted in panels (a), (b) and 
(c), respectively. In all three panels the initial values are drawn with dash-dotted lines, at t = 2.1 Myr with dashed lines, at t = 3.4 Myr 
dotted lines and at t = 4.2 Myr with solid lines. 
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Fig. 3. — Comparison of the temperatures corresponding to intermediate densities under the thermal-equihbrium (dotted hue) and isobaric 
(dashed line) regimes. The leftover accreting gas at the end of the condensation process flows along the nearly isobaric regime. 
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Fig. 4. — Pixel by pixel plot of log(P/fc)-log(p) for run DENS at different times. For reference, the equilibrium log(P/fc) curve (see Fig. 1) 
is also plotted. Time is shown at each panel in units of 1 Myr. Note that two of the diagrams were taken at the same time than the snapshots 
of Fig. 2. 
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Fig. 5. — Same as Fig. 2 but for run DEN15. In panel (c) the curve labeled with Cs represents the adiabatic sound speed at 6.2 Myr. 
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Fig. 6.— 



Same as Fig. 4 but for run DEN15. Note that three of the diagrams were taJjen at the same time than the snapshots of Fig. 5. 
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Fig. 7. — Time development of run DEN75. The density, pressure and velocity profiles at different times are plotted in panels (a) and (d), 
(b) and (e), and (c) and (f), respectively. The arrow in panel (a) marks the location of a peak in density originating from the small-scale 
components of the initial gaussian fluctuation (see text). 
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Fig. 8. — Same as Fig. 4 but for run DEN75. Note that the diagrams were taken at the same time than the snapshots of Fig. 7. 
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Fk;. 9. — Profiles of tlic density, sound speed (divided Ijy 4), velocity, and cooling time profiles, between and —75 pc at t = 27 Myr for 
run DEN75 (panel a). Panel (b) shows a zoom of the central parsec of this run. The density {solid line) is in cm~'^, the velocity {dashed line) 
and the adiabatic sound speed (triple-dot dashed line) are in km/s, P/k {dotted line) is in K cm~^, and the temperature {dot-dashed line) is 
in K. Note that the simulation is symmetric with respect the origin. 
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Fig. 10. — Comparison of density (panel a), pressure (panel b) and velocity (panel c) for run DEN75 with different resolutions. In panel 
(b) the pressure for the low-resolution simulations has been shifted 0.05 to make the plot readable. 
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Fig. 11. — Density {solid line) and velocity (divided by 3) {dot-dashed line) profiles at different times for runs SP20a {left panels) and SP4a 
{right panels). In the upper panels the initial density configuration is also plotted {dotted line). Clouds undergo an overshoot of density in 
the crushing stage (not shown). Most of the fully-condensed clouds present densities between 40 and 80 cm"^ but can be larger when two 
clouds collide. 
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Fig. 12.— 



Same as Figure 11 but for runs SP20b {left panels) and SP4b {right panels). 
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Fig. 13. — Density {solid line) and velocity {dot-dashed line) profiles for runs SP20c,d,e {left panels) and SP4c,d,e {right panels). Note that 
for spectra with femax/^min = 4 the evolution proceeds more slowly. 



Table 1 

Summary of runs with a density spectrum. 



Run 


S 


I [pc] 


^max / ^min 


SP20a 


3/2 


250 


20 


SP4a 


3/2 


250 


4 


SP20b 


3/2 


25 


20 


SP4b 


3/2 


25 


4 


SP20c 


1/2 


250 


20 


SP4c 


1/2 


250 


4 


SP20d 


1/2 


25 


20 


SP4d 


1/2 


25 


4 


SP20e 


1/2 


2.5 


20 


SP4e 


1/2 


2.5 


4 



Table 2 

Summary of runs with forcing. 



Run 


I [pc] 


mesh points 


/o 


(e) 


Mrms 


'^cond 


1 


3 


1200 


0.01 


0.985 


0.03 


11.9,12.1 


2 


3 


1200 


0.03 


0.990 


0.07 


8.4,10.85 


3 


3 


1200 


0.09 


1.095 


0.18 


6.15,7.8 


4 


3 


1200 


0.27 


1.880 


0.30 


3.95,4.6,6.40 


5 


3 


1200 


0.36 


2.450 


0.30 


> 16.0 


6 


3 


2400 


0.36 


2.530 


0.30 


> 12.0 


7 


3 


2400 


0.48 


2.610 


0.30 


> 14.0 


8 


3 


600 


0.72 


3.070 


0.40 


> 26.0 


9 


3 


3600 


1.20 


3.300 


0.65 


oo 


10 


10 


1200 


0.01 


0.980 


0.027 


10.80,11.20 


11 


10 


2400 


0.04 


0.980 


0.10 


7.30, 7.40 


12 


10 


7200 


0.04 


0.980 


0.10 


7.45,8.60 


13 


10 


1200 


0.10 


0.995 


0.13 


5.1,5.4,7.0 


14 


10 


7200 


0.18 


1.035 


0.25 


3.75,3.90 


15 


10 


1200 


0.24 


1.135 


0.30 


4.4,4.6 


16 


10 


1200 


0.72 


1.945 


0.60 


1.55,1.70 


17 


100 


9200 


0.09 


0.940 


0.30 


1.78,1.90 



